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ABSTRACT 

f—i Disk accretion onto weakly magnetized astrophysical objects often proceeds via a 

^ boundary layer that forms near the object's surface, in which the rotation speed of the 

^ accreted gas changes rapidly. Here we study the initial stages of formation for such a 

D boundary layer around a white dwarf or a young star by examining the hydrodynamical 

shear instabilities that may initiate mixing and momentum transport between the two 
fluids of different densities moving supersonically with respect to each other. We find 
that an initially laminar boundary layer is unstable to two different kinds of instabilities. 
p^ ' One is an instability of a supersonic vortex sheet (implying a discontinuous initial profile 

^ of the angular speed of the gas) in the presence of gravity, which we find to have a growth 

^ (-H rate of order (but less than) the orbital frequency. The other is a sonic instability 

Oh of a finite width, supersonic shear layer, which is similar to the Papaloizou-Pringle 

O instability. It has a growth rate proportional to the shear inside the transition layer, 

^ which is of order the orbital frequency times the ratio of stellar radius to the boundary 

Ci layer thickness. For a boundary layer that is thin compared to the radius of the star, 

the shear rate is much larger than the orbital frequency. Thus, we conclude that sonic 
instabilities play a dominant role in the initial stages of nonmagnetic boundary layer 
^ formation and give rise to very fast mixing between disk gas and stellar fiuid in the 

Q supersonic regime. 

cn 

^ Subject headings: accretion, accretion disks - hydrodynamics — waves - instabilities 



1. Introduction 

Accretion onto astrophysical objects possessing a material surface (as opposed to accretion 
onto black holes) always involves a non-trivial interaction of the incoming gas with the object's 
outer layers. Examples of such objects include white dwarfs in cataclysmic variables (CVs), young 
stars gaining material from a protoplanetary disk, and accreting neutron stars. In all these systems, 
one must understand how the incoming gas shares its angular momentum with the accreting object 
and how it mixes with the previously accreted material. 
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If the magnetic field of the central object is strong enough, it can disrupt the disk at some 



distance above the surface. Inside this region, accretion proceeds along field lines (Ghosh & Lamb 



1978 Koldoba et al. 2002). The critical value of the magnetic field at the surface of the central 



object for magnetic disruption is given by 
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(1) 



Here, we have used parameters for the mass and radius that are typical for CVs in outburst 



(Bergeron et al. 1992; Hachisu 1986 Patterson 1984). The parameter /3 is a dimensionless factor 



of order unity that depends on the model adopted for the disruption of the disk by the stellar 



magnetic field (Ghosh & Lamb 1978). 



For B^<B^, 

crit) the accretion disk will extend all the way to the surface of the star, for which 
there is good observational evidence in neutron star low mass X-ray binaries ( Gilfanov et al.|[2003 ) 



and in dwarf nova systems in outburst (Wheatley et al. 2003). Since the rotation rate of the star, 
f]*, is slower than the Keplerian rotation rate at the stellar surface, ^k{R*)^ a boundary layer 
(BL) will exist inside of which dO,/dR > and the rotation profile of the star attaches smoothly 
to that of the disk. If the accretion rate is high enough, it is also possible for material to spread 



meridionally to high latitudes forming a spreading layer (SL) (Inogamov Sz Sunyaev 1999; Piro &; 



Bildsten 2004). For both a BL and a SL, the intense energy release in a localized region near the 
stellar surface leads to easily observable signatures such as hard spectral components, variability 
of emission, and so on. 

A considerable amount of effort has been previously devoted to understanding the structure of 
well-developed, steady-state BLs, which have been evolving for long enough to establish a smooth 
rotation profile in their interior. An outstanding question in the study of steady-state BLs is iden- 
tifying the mechanism of angular momentum transport in the layer. Several potential mechanisms 



Kippenhahn & Thomas 


1978), 


er-Spruit dynamo ( 


Spruit||1999 



Piro & Bildsten 2007). Despite this effort, no clear answer exists at present regarding the nature 



of the angular momentum transport and mixing in well-developed BLs. 



1.1. Initiation of Boundary Layers 

An interesting aspect of the BL problem that has received less attention in the past is the issue 
of BL initiation, i.e. the initial stage of BL formation which must occur when the accreted material 
first touches the surface of the star. In this case the physical setup is going to be quite different 
from the steady-state BL primarily because of the much larger velocity shear, justifying the study 
of BL initiation in the limit of an essentially discontinuous rotation profile at the interface between 
the star and accreting material. 
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Even though the initiation stage is just a transient phase in the BL evolution, one still needs 
to understand it to get a full picture of the BL phenomenon. Also, one need not think that the BL 
initiation is a unique event for every accreting object — it may, in fact, be repetitive. The most 
common situation in which this BL initiation is recurrent involves the dramatic increase of the mass 
accretion rate (due to some sort of outburst triggered by an instability) through the disk which is 
magnetically disrupted outside the star under normal circumstances. According to Equation ([T]), 
a sudden increase of M by several orders of magnitude as typical for some accreting systems (e.g. 
dwarf novae, FU Ori outbursts, etc.) can rapidly compress the stellar magnetosphere to the point 
at which disk material starts touching the stellar surface and a BL starts to form. 



There is good observational evidence for this kind of recurrent behavior. For instance, Livio 



Sz Pringle ( 1992 ) have argued that the observed lag in the rise time of the UV emission relative to 
the optical in a CV system transitioning to outburst can be explained if the disk is magnetically 
disrupted in quiescence but not in outburst. Thus, the hottest, innermost part of the disk is 
evacuated in quiescence, and UV radiation is not emitted immediately in the transition to outburst, 
since the empty region requires time to become filled. 

FU Ori stars are another type of system in which the disk can extend all the way to the stellar 
surface, due to the high accretion rate - M ~ lO~^M0yr~^ ( Kenyon et al.|[r988 ). In these systems. 



the BL can puff up to become of order the stellar radius and the distinction between the boundary 



layer, the disk, and the star becomes blurred (Popham et al. 1993) 



The goal of our present work is to make a first step towards understanding the formation of 
the BL. As in the case of a well-developed BL, the major issue for this initial stage lies in details of 
the angular momentum transport and mixing. However, because of the enormous shear present at 
the star-disk interface at this stage, it is highly likely that purely hydrodynamical shear instabilities 
would dominate the transport of mass and momentum rather than anything else. 

For that reason, in this work we primarily focus on exploring the properties of various shear- 
driven instabilities under the conditions typical at the BL formation stage. First, we seek to identify 
a particular variety of the shear instability that most efficiently initiates mixing between the two 
fluids, i.e. has the fastest growth rate. Second, we examine the conditions needed for its operation, 
such as the density contrast between fluids, initial velocity profile, etc. 

There are several physical ingredients which can potentially affect operation of shear insta- 
bilities: stratification, rotation, magnetic fields, turbulence, radiative transfer, and the supersonic 
nature of the flow. The latter aspect of the problem is very important and is inevitable during 
the BL initiation phase, when disk material rotating at the Keplerian speed comes into contact 
with the more slowly spinning stellar surface. The differential azimuthal velocity between the two 
interacting flows is then bound to be a significant fraction of the Keplerian velocity at R^: and 
should highly exceed the sound speed both in the disk and in the outer layers of the star. 
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1.2. Shear Instabilities in Compressible Fluids 



The study of shear instabilities in compressible fluids is of fundamental physical significance and 



has received much attention in the past. Landau (1944), Hatanaka (1949), Pai (1954), Miles (1958), 



and Gerwin ( 1968 ) have all studied the problem in the vortex sheet approximation, where one half 



plane of compressible fluid moves at constant velocity over another. When the two fluids move at 
a low Mach number relative to each other, one finds that infinitesimal perturbations are governed 
by the classical Kelvin-Helmholtz (KH) dispersion relation for two incompressible fluids. However, 



Miles (1958) showed that above a critical Mach number, the vortex sheet becomes marginally 



stable to two-dimensional disturbances along the direction of the flow. This is a surprising result 
considering that one might have expected increasing the shear to lead to an increased growth rate 
for the instability rather than stabilization. Understanding the implications of this result for the 
momentum and mass transport in the astrophysical BLs is one of the goals of our study. 



Later, Blumen et al. (1975), Ray (1982), Choudhury & Lovelace (1984) and Glatzel (1988) 



studied the stability of a fluid with a continuous and monotonically varying velocity proflle. They 
found that unlike the vortex sheet, continuous velocity proflles with a supersonic velocity difference 



across them were unstable even at high Mach number. Glatzel (1988) showed that the instability 



was similar to the Papaloizou-Pringle instability which operates in hydrodynamical disks with radial 
boundaries ( Papaloizou &: Pringle|[T"984 ) . Thus, for high Mach number flow in a compressible fluid, 
a finite thickness velocity profile exhibits fundamentally different behavior from the vortex sheet, 
since it can support unstable modes. Moreover, the growth rate of the unstable modes scales as 
(X 1/L, where L is the thickness of the shear layer. Thus, the thinner the shear layer, the faster the 
instability operates! This is exactly the opposite of what one might have expected from the vortex 
sheet stability criterion and we will provide an explanation for this in ^ 

Previous studies of possible instabilities inside astrophysical BLs have been primarily concerned 



with the sub-sonic regime when the flow can be considered as almost incompressible (Kippenhahn 



fc Thomas|1978HFujimoto|1988D . This regime may indeed apply in well-developed BLs with smooth 
shear proflles, even though at the present level of knowledge one can hardly exclude the possibility 
of the existence of localized regions in steady-state BLs where compressibility is still important. 
Our primary goal here is to extend these studies into the regime of highly compressible, supersonic 
flows and to explore the implications for astrophysical objects. In the course of studying supersonic 
shear instabilities, we will sometimes also incorporate stratification in our calculations and explore 
the role of rotation. 

The paper is organized as follows. In ^we present our governing set of equations and describe 
the formalism we use to study the stability problem. Then in ^we study the case of the compress- 
ible vortex sheet with no gravity and show that our formalism reproduces the dispersion relation 



obtained by previous authors (Landau 1944 Hatanaka 1949 Pai 1954 Miles 1958 Gerwin 1968). 
In we introduce gravity as a small parameter and perform a perturbation expansion to study 
what influence this has on the stability of the vortex sheet. Finally, in ^we study the stability of 
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a finite width shear layer with a linear velocity profile. We derive new dispersion relations for this 



case and test the growth rate of the fastest growing mode using the Godunov code Athena (Stone 
et al.||2008p. 



2. Formalism 

Here we discuss the formalism that we use for analytic calculations throughout the paper. We 
adopt cylindrical coordinates {zu, 9, z) and for simplicity assume that all equilibrium quantities are 
independent of z. Thus, for our purposes gravity is given by —g{vu)'d7 (giw) > 0), and we ignore 
the vertical stratification of the disk. Such a setup does not allow for the baroclinic instability, 
which would require dVl/dz ^ 0, but it does contain the necessary ingredients for studying shear 
instabilities. Given our assumptions, the equation of hydrostatic equilibrium reads 

1 dP 

g{w) = g{w)-wVL^{w) = — — , (2) 

p dw 

where g{vj) is the effective gravity. The equilibrium density, pressure, and sound speed are given 
by pi'ccj), -P(ro), and s{w) respectively. 

Denoting the {w^O^z) velocities by {u^v^w), using v = Vlw, and assuming adiabaticity, the 
Euler equations in cylindrical coordinates are 

dp Id, ,ld,,d,, 

ot w ow w oO oz 

Du IdP 
Dt p ovj w 

Dv 1 dP uv 

'Dt ~ ^~pw^ ~ ^ ^ ' 



Dw _ IdP 
~Dt ~ ~~plh 



(6) 



D{Pp~^) 



^ ^ (7) 

Dt ^ ' 

D d d V d d , . 

Dt dt dw w dO dz 

On top of the zeroth order equilibrium state, we consider infinitesimal two-dimensional pertur- 
bations of the form df{w) exp[i{m9 + k^z — ut)], where the 6 denotes an Eulerian perturbation. 
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Starting from Equations ([3])-([7|), the linearized first order equations are given by 



iu5p 
ioj5u 

iCoSv 

icoSw 

iuSp 



w 
1 



-{vjpbu)' 



imp 
w 



6v + ikzp5w 



6P' + ^6p 
P P 



6P + 2B6u 



im 
vjp 

P 

-Clp6u + ^6P, 



(9) 
(10) 
(11) 
(12) 
(13) 



wliere w = w — mil is the phase speed in the locally corotating frame, B = {Vt + d{'UjQ) /dw) /2 
is Oort's B constant, Cl = ^^^dlnP/dw — dlnp/dw is the Ledoux discriminant, and the primes 
denote differentiation with respect to w. Note that in Equations ([9])-(13), we have made Cowling's 



approximation by ignoring the first order perturbation to the gravitational potential (Cowling 



1941). This means that we are taking the self-gravity of the BL to be negligible. 



We can reduce the above set of linear equations to a pair of equations in 6P and 6u: 

2Bm\ 



ki 



m 



5P 



P ClCo + 



UJ 



5u H (zupdu)' 

w 



ip{uj — gCi — K )5u 
Here we have introduced the epicyclic frequency 



OJOP H ttOP 



2Q.m 



-6 P. 



(14) 
(15) 



ABn. 



One generally expects the BL width 6bl ^ohe small compared to the stellar radius, 5bl ^ R*- 
Under this assumption, we show in Appendix |A] that Equations (14)-(15) can be simplified to the 
following set of equations: 



OJ 



m 



6P 



ip{uP' - gCL)Su 



u6P' + ^5P 



Su + ijj{p5u)' 



(16) 
(17) 



The form of Equations ( 16 ) and ( 17 ) is the same as for a plane-parallel stratified shear flow ( Alexakis 



et al. 2002). This suggests that the effects of Coriolis force and curvature are unimportant in a 



radially thin BL if shear instabilities provide the turbulence, and leads us to redefine the cylindrical 
{w, 6, z) coordinate system into a Cartesian (x, y, z) system. The perturbed quantities now have 
the form (5/(x) exp[i(fcj^?/ + kzZ — wt)], where ky = m/R^,. From this it immediately follows that 
u) = OJ — kyVy{x), where Vy{x) = i?*r2(x) is the velocity profile. From here on, we will ignore the 
rotation terms, and treat the flow in the BL as a plane parallel shear flow. 



From Squire's theorem for plane parallel shear flows (Fejer & Miles 1963), any three-dimensional 
perturbation is mathematically equivalent to a two-dimensional perturbation upon making the 
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transformations ky ^ ky/ cos 9, A;^ — )• 0, and Vy — )• VyCosO, where cos 6* = ky/ ^k'^ + k1. Thus, it 
is sufficient to only consider two-dimensional perturbations {k^ = 0) in the stability problem, and 
for the rest of this paper we take perturbations in the form 



5f{x) exp[i(fcj/y - ut)]. 



(18) 



Assuming two-dimensional perturbations. Equations (16) and (17) become 

,2 



kl 1 5P 



p {ClO + kyVy) 5u + 0{p5uy 



ip{u? - gCL)5u 



udP' + ^6P, 



(19) 
(20) 



where the primes now denote differentiation with respect to x. Equations (20) and (19) can be 



used to obtain the generalized Rayleigh equation (Appendix B.l), which using a notation similar 
to 



Alexakis et al. ( 2002 ) reads 



b" + kl 



_kg -\- kq 

qp — = — - 




0. 



Here, 



-6u^/p/kx, 



(21) 



(22) 



is a modified stream function, and for simplicity we have dropped the bar over so now g denotes 
the effective gravity. We have also defined the quantities: 



W = kyW^p/ik^ 

W = Vy — c where c = oj/k is the phase speed. 

ky{W'^ /s"^ — 1) is the square of the x-component of the wavevector 
in the absence of shear or stratification. 
p = pf"^ , where p is the density. 



krf 



f = exp 



kgiOdC 



(23) 
(24) 
(25) 
(26) 
(27) 

(28) 

(29) 
(30) 



kg = g / \s a. measure of the inverse of the local scale height. 
ks = p I p is the inverse stratification length scale. 

Note that the Ledoux discriminant is given by Cl = —(kg^kg). 

The boundary conditions on (50 that must be satisfied at a discontinuity in the density or the 
velocity are that the upper and lower fluids stay in contact and that the pressure perturbation is 



continuous across the interface. We show explicitly in Appendix B.2 that these conditions can be 
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formulated as 

^ = ^ (31) 

W+ W- 

W+5<t^'^-W'^5<t^+-^-±h^ = W^S<p' -W'6c^^-^-^h^. (32) 

W+ W- 

The +/— signs denote evaluation of a quantity in the upper/lower fluid at the location of the 
interface. 



3. The Vortex Sheet without Gravity 



Equations (21)-(32) are fully general and apply to an arbitrary velocity profile. In particular, 
we can assume that the velocity varies discontinuously at some radius x = 0: 



Vyix) 



Vy, X>0 
-Vy, X <0, 



(33) 



where is a constant. This is known as a vortex sheet approximation. It represents the simplest 
possible description of the velocity variation between the two limiting values by essentially ignoring 
the details of the transition. Subsequently in ^ we explore a more realistic model of the velocity 
variation, in which the transition occurs in a region of finite radial width. We point out that the 



vortex sheet approximation is valid for kydsL ^ 1- We show in Appendix C.3 that in this limit the 



dispersion relation for a finite width layer of constant shear reduces to the vortex sheet dispersion 
relation. 

We assume in this section that 5 = 0, so there is no gravity, and that p and s are constant 
above and below the interface, but can be discontinuous across it. This case has already been 



considered by other authors in the past including Landau (1944), Hatanaka (1949), Pai (1954), 



Miles (1958), and Gerwin (1968). However, the results for the case without gravity will often be 
referenced later in the paper and serve as a verification of the formalism we developed in ^21 



Using the velocity profile (33), the generalized Rayleigh equation (21) becomes 



(34) 
(35) 



where just as in 1§ the +/- signs denote the upper/lower fluids respectively. Since kx,± is constant 
for Vy{x) given by Equation (33), we have that 



oc 



(36) 



In general, kx^± is complex, and the sign is determined by applying the appropriate boundary 
conditions. We discuss the boundary conditions shortly, which will also make clear the reason 
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for the negative sign in the exponential of Equation (36). Plugging the expressions for 5(l)± into 



Equation (32) gives 



W-kr-C 



(37) 



where we have used 5(f)'^ = ikx,±S(j)±, p± = p±, and W!^ = 0. Using Equation (31) to substitute for 
S(j)- in terms of (5(/>+ , and substituting for W in terms of W and p we have 



kx,- 



k 



(38) 



Introducing the density ratio e = p^j p-^ the Mach number in the upper fluid M = V^/s+, and the 
phase speed normalized by the sound speed in the upper fluid 99 = c/s+, we have 



e(M - ip)\ {M + (^)2£± - 1 = (M + (^)V(^- 95)2 _ 1, 



(39) 



By the definition of the sound speed, = 7-P/p, the condition of pressure balance everywhere 
throughout the flow requires that j^^p-^s\, = ^Z^ P-S^_, where 7+ and 7_ are the adiabatic indices 
above and below the interface. Assuming 7-|_ =7-, we have from pressure balance that {s-/ Sj^Y — 
p+/ P- = e. Thus, the dispersion relation becomes: 



e(M - ^fy/{M + ^Ye-^ - 1 = (M + ipfyJ{M-ipf 



1. 



(40) 



Miles ( 1958 ) has studied Equation \Am and found that the stability criterion, i.e. that ip is purely 



real, is given by 



M > Merit = ^(1 + 6^/^/'. 



(41) 



This shows the surprising result that infinitesimal disturbances are stabilized at high Mach number. 



However, Fejer & Miles ( 1963 ) pointed out that due to Squire's theorem, it is always possible to 



choose an angle 9 for the wavevector with respect to the flow velocity such that the projected Mach 
number M cos6 is smaller than Merit- Thus, even at high Mach number, the vortex sheet without 
gravity is still unstable to three dimensional disturbances, which are almost perpendicular to the 
direction of the flow; these unstable oblique modes resemble classical KH modes. However, in the 
astrophysical context, a thin disk has a scale height s/Q ^ i?*, and the wavelength of the oblique 
modes will be small relative to the disk scale height only for very small wavelengths X <^ s/0,. If 
the BL itself has a thickness 5bl ^ s/fl, the vortex sheet approximation for the oblique modes is 
invalid, since either the modes don't fit into a disk scale height, or the condition A ^ 6bl is not 
satisfied. 
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3.1. Solutions for M > 1 



Since we are interested in the high Mach number hmit for the initiation of the BL, we now 



find analytical solutions to the dispersion relation (|40|) for M ^ Merit, where M^it was defined in 

(42) 



Equation (41). Squaring Equation (40), one obtains a sixth order polynomial in ip: 

(g2 + ^fe){M - iff - (1 - (M - vf){M + iff 



0. 



This polynomial has two easy to find analytic factors (|Miles 1958|) 

.1/2 \ 



-M 



1 



-M 



1 + 



(43) 



We now take the limit M ^> Merit, in which case to terms of order 0{M~'^), the other four solutions 



to the polynomial in Equation (42) are 



M + 1 + 



2(2M + 1)2^ 



M- 1 



M + + 



2(2M - 1)2' 

6^/2 



2(2M-ei/2)2' 



-M 



.1/2 



.1/2 



2(2M + ei/2)2 



• (44) 



It is easy to see that the six roots of the polynomial (42) correspond to sound waves. Starting 



from Equation (25) and rearranging terms we obtain 



(45) 



As long as kx is real, this is the dispersion relation for a sound wave, since W"^ is the square of 
the phase velocity in the frame comoving with the fluid. Plugging in the six roots from Equations 
(43) and (44) into Equation (25), it is straightforward to verify that kx± are indeed real for all 



of them. Each of the four roots in Equation (44) has a further simple interpretation. The first 
two correspond to sound waves that propagate almost parallel to the interface in the +y and —y 
directions in the upper fluid, whereas the second two correspond to sound waves that propagate 
almost parallel to the interface in the +y and —y directions in the lower fluid. The two roots in 



Equation (43) are more difficult to interpret, but the first of these corresponds to a standing wave 
when the two fluids have equal density (i.e. e = 1). 

Since each of the six roots for M ^ Merit yields a real value for kx^±, the solutions do not 
damp away from the interface, and we need to apply radiation boundary conditions at x = ±00. 
The proper procedure is to demand that all waves are outgoing in each fluid in a frame which 



is subsonic with respect to the fluid (Miles 1957b), and it is convenient to work in the comoving 



frame of each fluid. The dimensionless phase velocity in the frame comoving with the upper fluid 
is pi+fiF = P — M, and the analogous expression for the lower fluid is p>-,cF = f + M. 
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In order to satisfy Equation (38), kx^-\- and k^^- must have the same sign, so we must have 
either 6<f>± oc e*(l'^^'±l^'~'^*) or 6(f)± oc e'^i-\kx,±\x~uit) ^ Moreover, since kx,+ and kx,- have the same 
sign, it is clear that (p+,cF and (p-,cF must have the opposite sign to yield outgoing waves in the 
comoving frames of each of the two fluids. Only three of the six roots found above satisfy this 
condition: 



2(2M-ei/2)2 



-M 



1 



=1/2 



1 + el/2 



M - 1 



2(2M - 1)2' 



(46) 

(47) 
(48) 



We will refer to these three roots as the lower, middle, and upper branches, respectively. Further- 
more, it is straightforward to check that the solutions which yield outgoing waves in both the upper 
and lower fluids have kx ± are positive given our definition d36^ . 



3.2. Dispersion Relation in the General Case 



We now relax the assumption of M ^> 1 and numerically solve the dispersion relation ( 40 ) at 



arbitrary Mach number for e = 1 (Figure [T^) and e = .01 (Figure [Tja) . At the critical Mach number 



given by Equation (41 ), the upper and lower branches merge together in the real plane and bifurcate 
in the complex plate. These bifurcated solutions turn into the two incompressible KH modes for 
M <^ Merit- Unlike, the lower and upper branches, the middle branch has no incompressible analog 
and ceases to be a viable physical solution below a critical Mach number 

2M^ = l + ei (49) 



The reason for this is that below M = Mm, kx,± switches from real to imaginary and the boundary 
conditions for the middle branch at rr = ±00 can no longer be satisfied. 



4. The Isothermal Vortex Sheet with Gravity 

We now go beyond the simplifying assumption of no gravity used in ^ and consider the case 
where g is non-zero and constant. We shall shortly assume an isothermal equation of state, but for 
now we consider the more general polytropic equation of state of the form P = K±p'^ in each of 
the two fluids. For any equation of state of this form, we have 

dPdp , - 

^, = 4 = -^ = --^ = --^^> m 

'J dr 7 
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Fig. 1. — The dispersion relation (p/M as a function M for e = 1 and e = .01. The soHd and 
dashed curves give the real and imaginary components respectively. The blue, black, and red 
curves correspond to the lower, middle, and upper branches respectively. The red and blue curves 
have been slightly offset vertically so they do not overlap. 



which means kg + kg = (1 — j/n)kg. If n = 7, corresponding to the case of an adiabatic atmosphere. 



the second term in parentheses in Equation (21 ) drops out, but then the third term becomes difficult 
to treat analytically. However, for an isothermal atmosphere (n = 1) the sound speed is constant. 



and both the second and third terms in Equation (21 ) reduce to a tractable form. Since we are only 



interested in the gross, qualitative properties of the flow, we assume both fluids are isothermally 
stratified, since this assumption significantly simplifies the analytical treatment. 



For n = 1, the second term in Equation (21) becomes 

'^9 



^kg -\- kq 
qp — - — 
^'^ 11/2 



klW^ 



(51) 



We now again assume a vortex sheet velocity profile (Equation (33)), in which case for n = 1, the 

■2-7, 



third term in Equation (21) is given by 

W" (Vp)" 



w 



(52) 



Putting Equations (51) and (52) into Equation (21), we have 

b" + kl5^ = 0, 



l-(7-l) 



kg\ /^\2 



kyj \WJ 



(53) 



The perturbation is then given as 



t>± (X e 



(54) 



We now comment on the terms present in Equation (53). In the absence of gravity, kg = 0, 
and we simply have kx = k^- If kg ^ 0, then the second term on the right hand side of Equation 



(53) can be written in a more famihar form as 

2 



(7-1) 



(- 



N 



UJ — kyVy 



where 



(7-1)5' 



(55) 



(56) 



is the Brunt- Vaisala frequency. Thus, the second term on the right hand side of Equation ( |53| ) 
provides a lower frequency cutoff for sound waves at the Brunt- Vaisala frequency. The last term 



in Equation (53) arises because we have not made the short wavelength approximation k^ ^ kg. 



Now, we use the boundary conditions (31) and (32) to determine the dispersion relation. 



Substituting 6(j)'j. = —ikx^±54>±, and VFj_ = (2 — ^)kg^±W±/2 into the condition (32), we have 



7 g 



+ 



7 9 



9P~ 



(57) 



Next, using the second boundary condition (31) to substitute for 5(f)^ in terms of 



the fact that p± = p± at the interface, we have 



7 9 



Wi 



9P+ 



2 

We now introduce the dimensionless gravity parameter 

9 



+ 



7 9 



Wt - gp. 



G 



k s2 ' 



and using 
(58) 

(59) 



which is closely related to the ratio of the wavelength to the pressure scale height, hg. Thus, G ~ 1 
when kyhs^+ ~ 1, and Ge^^ ~ 1 when kyhg^^ ~ 1. Using e = p+/p- = s^_/s\, and performing 
some algebra, we have 



. kr,- 



kii 



+ 



7 



G 



kn 



e + G(l-e) 



. krr 



kij 



^Ge-' 



ko. 



kx. 



2 fW^ 



Setting G = it is clear that we recover the vortex sheet dispersion relation in the absence of 
gravity (Equation (p8|)). 



We now check Equation ( 60 ) by showing that it reproduces the well-known incompressible KH 



dispersion relations in the limit M <^ Merit, before going on to treat the case of the supersonic 
vortex sheet with gravity. 
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4.1. Highly Subsonic Vortex Sheet with Gravity 



We assume that V / s± <C 1, which immediately implies M <C Merit, and we also assume G <C 1 
and Ge^^ <^ 1, which means that the wavelength of the perturbation is much smaller than the scale 
height in both the upper and lower fluids [kg^±/ky <^ 1). Next, we eliminate sound wave modes by 
assuming ip 1 (phase velocity much lower than sound velocity in upper fluid) and (/9e~^/^ ^ 1 



(phase velocity much lower than sound velocity in lower fluid). According to Equation (53), this 
means kx ^iky^ and given our definition for 64>± in Equation (153), we must chose the minus 
sign in the upper fluid and the plus sign in the lower fluid to give vanishing solutions at x = ±00. 



Equation ( 60 ) then yields 



UJ 

kii 




(7 



(M 



-^ + ^(1 




(7-l)GV 



(61) 



Next, we note that (7 - 1)0"^ /{M - ipf 



kyVyf and that (7 - l)G'^e~^/{M + <ff 



N'^/{uj + kyVy)'^, where is the Brunt- Vaisala frequency for an isothermal atmosphere and was 



given in Equation (56). We expect \uj ± kyVy\ > Vgk, which is the characteristic frequency of 

\/9kg,±, and we have already assumed kg^±/ky <C 1, it follows 



surface gravity waves. Since A^it 
that N± ^ |a; =F kyVy\. Consequently, Equation (61) reduces to 

2 



— e 



V, 



y 



(62) 

%y / f^y V^y 
which is the well known KH dispersion relation for an incompressible fluid in the presence of gravity. 



4.2. The Weak Gravity Limit at High Mach Number 



As mentioned before in [ 3.1 Miles ( 1958 ) has demonstrated that in the absence of gravity, the 



vortex sheet is stable above a critical Mach number given by Equation (41). We now address the 
question of whether the system still remains stable when G 7^ 0. 



To answer this question, we will use our general dispersion relation (60) in which we will 



additionally assume M ^ Merit, since this assumption significantly simplifies the algebra. Since 
we have already obtained solutions for the case G = and Ai S> Merit ( ^3.1[ ), we proceed by 
considering gravity as a perturbation. We consider wavelengths that are much smaller than the 
pressure scale height and ask what happens in the limit G — )■ 0. In our analysis, we will assume 
that the density ratio e < 1, so that the system is stable to the Rayleigh- Taylor instability. 

Because the introduction of gravity makes the upper and lower fiuids stratified, some care 



should be taken when determining which solutions are physical and which are not. Miles ( 1957b ) has 
shown that for G = 0, the three supersonic KH modes can be understood in terms of sound waves 
emitted from the interface between the two fluids. This interpretation is useful as well for the case 
with gravity and leads to a couple of insights. First, the amplitude of sound waves is not constant 
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as they propagate through a stratified medium. Rather, to conserve energy, waves propagating 
upward (to lower densities) must increase in ampHtude, and those propagating downward (to higher 
densities) must decrease in ampHtude. Second, since the sound waves are emitted from the interface, 
if CO has an imaginary component, then the amplitude of the emitted waves changes with time. Thus, 
if ^[u}] > 0, and there is an instability, then the amplitude of the waves will decay with distance 
from the interface, since the waves emitted in the past had lower amplitude. Conversely if '^[oj] < 
and the perturbation is damping in time, then the amplitude of the emitted waves will increase 
with distance from the interface, since the waves emitted in the past had a larger amplitude. Both 
of these effects mean that the amplitude of the waves can blow up as we move away from the 
interface. Thus, we take as physical those solutions which yield outgoing waves (away from the 
interface) in both the upper and lower fluids, even if these solutions diverge as x — )• ±00. For the 
limit G — )• 0, this means that the physical solutions are still the lower, middle, and upper branches 
but now modified by the presence of a weak gravitational field. 



We note here that although the dispersion relation (60) is valid for all values of G and not just 



for small G, the simple picture of outgoing sound waves in the upper and lower fluids is only valid 
for G — )• 0. From a physical point of view, this can be attributed to the following fact. Sound waves 
(p-modes) traveling in a stratified atmosphere have a frequency cutoff at cj^ = N'^ below which 
propagation is not possible. Given our definitions of ky and kx, the frequency of a sound wave in the 
frame comoving with the fluid is ~ {ky + kDs^. Substituting A^^ ~ and using the definition 



of A^^ (Equation (56)) yields < (7 — 1) -"^(1 + {\kx\/ky)'^) for sound waves to propagate. If this 



condition is not fulfilled, then the picture of outgoing sound waves is invalid, and the boundary 
conditions need to be formulated in a different way, which is beyond the scope of the present work. 



4.2.1. The Lower Branch 



We begin by considering how the lower wave is modified in the limit G — t- and M ^ Merit- 
If G = exactly, then kx = kx, and in the limit G — )• we have kx = A;^, (l + ©(G^)). Keeping 



terms only to first order in G, Equation (60) becomes 

kfr 



kit 



e + G(l-e) 



kx,+ J V s+ 

Writing this out explicitly in terms of M and (p gives 

i{M-ipYe i{M + ipf 



kx- y V s+ 

{M-pfe 
(M + (/9)2e-i - 1 ~ (M -iff -I 



.(63) 



.(64) 



Next, we assume that gravity only weakly affects the dispersion relation and make the pertur- 
bative expansion 



(65) 



- 16 - 



where ipo is the solution for G = and M ^ Merit- 



For the lower branch, we can use Equation (46) for ipQ, and in the limit M ^ Merit we have 

-M + 

a/2 



ky 



(66) 
(67) 



2^-1 



1 



1 



2M 



a/2 ■ 



Defining 



Ml = 2M 



.1/2 



Equation (64) becomes 



^{Mi - v9i)2 - 1 V(e'/2 + (^i) 



2^-1 



2-7 

+ G(l-6) = — 



(ei/2 + ^i) 



2,-1 



(68) 

(69) 
{Mi-^ife 

{Ml - V9i)2 - 1 



It will turn out (Equation (|72|)) that (/?! is proportional to G, so working to first order in G is 
equivalent to working to first order in ipi . Equation ( 70 ) then simplifies to 

iMie{l + 2e~i/2^_^ 



i{Mi-ipi)e 



1 + Mfe-y^^i 



G{1 



2-7 



G Mf 



1 + 2e-i/2 



1 + 2Mfe 



2f-V2 



(71) 



Assuming \Mfe '^^'^ipi\ <^ 1, and keeping terms to leading order in M/, we can solve for ipi in terms 
of G. 

2-7 G . 



^1,1 



2 Miei/2 



(72) 



We immediately see two things from Equation (72). First, ipi is purely imaginary, and second, 
if 7 < 2, ipi is negative and the perturbation damps, whereas if 7 > 2, ipi is positive and the 
perturbation grows. For realistic equations of state, 7 < 5/3, so the lower wave always damps. 



.(70) 



4-2.2. The Middle and Upper Branches 

We can find an approximate solution for ipi in the limit G — )• and M ^ Merit for the middle 
and upper branches in much the same manner as for the lower branch. The first order correction 
for the middle branch is 

7^(1-6^2). 

^l,m = -^ ^2 (^3) 



and for the upper branch is 



^.2-7 G6V2 . 

'^^'"^ 2 2M-l" ^^'^ 



-17- 



Just as in the case of the lower branch, (pi is purely imaginary for both the middle and upper 
branches. We see that ipi is negative for the middle branch if e < 1. However, for the upper branch 
if 7 < 2, ifi is positive and if 7 > 2, ipi is negative. This is opposite from the lower branch meaning 
that for any 7 7^ 2 in the limit M ^ 1, one of the two branches is always unstable and the other 
one is damped. For a realistic equation of state, 7 < 5/3 so the upper branch is the unstable one, 
and the lower branch damps. 



4-2.3. Numerical Verification 



We verify Equations (72), ( |73[ ), and (74) numerically by solving the fully general dispersion 
relation (60) and comparing with our analytical estimate for the parameters M = 5, e = .5, 
and 7 = 5/3. We plot Q{(p) vs. G in for both our analytical solutions (solid lines) and the ones 
obtained numerically (dashed lines) in Figure [2j The analytical solution converges to the numerical 
one in the limit G — )• 0. 



Upper Wave 




(b) 



(c) 



Fig. 2. — Plots of Q{^p) vs. G for the lower (a), middle (b), and upper (c) branches respectively, 
using M = 5, e = .5, 7 = 5/3. The solid curves show the numerical solution obtained by solving 



Equation ( 60 ) and the dashed curves correspond to the approximate solutions from Equations ( 72 ) 



(73), and (74). 



5. Sonic Instabilities in a Finite Width Layer of Constant Shear 



The calculations presented in ^3|4 were performed for the velocity profile (33) featuring a 
discontinuity at some radius. We now consider a more complicated (and more realistic at later stages 
of the BL evolution) initial setup in which the velocity between the two fluids varies continuously 
within a narrow shear layer. Unlike the vortex sheet, the finite width shear layer without gravity is 
known to be unstable at high Mach number ( Glatzel|1988 Choudhury &: Lovelace|1984 Ray||1982 ), 
and in this case the growth rate of the instability scales inversely with the width of the shear layer 
and in proportion to the shear, S ~ VLkR*/ ^BL- 



In the following, we extend some of the findings of 



Glatzel 



( 1988 ) to study sonic instabilities 
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in a finite width shear layer without gravity and apply them to the problem of the BL initiation. 
The setup we consider has the velocity profile 



V{x) 



V, x> 6bl, 

Vx/5bl, -5bl <x< 6bl, 
-V, X < -6bl 



(75) 



and the density profile 



p{x) 




(76) 



Pressure equilibrium again requires that p+s^ = p-s^, which sets the sound speed everywhere in 
the fiow, and as before we have e = p+/ p- and M = V /s+. 



Although we only consider a linearly varying velocity profile in the shear layer, Ray (1982) 
and Choudhury & Lovelace ( 1984 ) have found that different shear profiles are qualitatively similar. 
Thus, we consider a constant shear to be representative of more general shear profiles. 



5.1. Dispersion Relation 



We now study the dispersion relation of the finite width shear layer. In applying the dispersion 
relation to the initiation of the BL, we are most interested in the growth rate of the fastest growing 
mode for M ^ 1. Glatzel (1988) has already obtained the dispersion relation for the case e = 1, 



and using his techniques, we derive the dispersion relation for the case of arbitrary e in Appendix 



[C| We also show in Appendix C.3 that the dispersion relation for a finite width shear layer reduces 
to the dispersion relation for a vortex sheet in the limit /cy — )■ at constant 5bl- 

In Figure [sj we plot the dispersion relation ( C20 ) as a function of wavenumber for the pa- 
rameters M = 5 and e = 1, e = .25, e = .1, e = .01, and e = 0. The e = case is equivalent to 
having a hard reflecting boundary at a; = —5bl- The modes depicted are the ones that converge 
to the upper and lower branches from ^ in the vortex sheet limit {kydBL ^ !)• For e < 1, the 
upper branch always has a larger growth rate than the lower branch, and both the upper and lower 
branches have the same growth rate for e = 1. In the e = 1 case, we can identify the upper and 
lower branches as the n± = decoupled modes in §5.4 of Glatzel (1988). In addition to the upper 



and lower waves, Glatzel (1988) has shown that there is an infinite spectrum of damped modes. 



but we do not consider these here, since we are interested in determining the growth rate of the 
fastest growing mode. 

We now revisit the seemingly paradoxical statement that the vortex sheet is stable to two 
dimensional disturbances along the flow direction above a critical Mach number; yet at the same 
time, the growth rate of the fastest growing mode scales with which implies that the thinner 
the shear layer, the faster the instability proceeds. The key to resolving this apparent controversy 
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is to consider a mode having ky5BL ^1- Its growth rate is diminished if we decrease 5bl while 
keeping ky constant for M > Merit and becomes vanishingly small if we take the limit 5bl 0. 
This can be seen from the dying left hand tail of the curve in Figure [Sja, and keeping ky constant, 
while decreasing 5bl we move leftward along the tail. This means becomes smaller, since 
ip = uj/kyS+ is directly proportional to u for constant ky even as we decrease 6bl- We next point 
out that the curve in Figure [3|: remains unchanged in shape or amplitude as we decrease 6bl, 
keeping kySBL constant. Consequently, the value of ky for which the maximum in '^[uj] occurs 
ky,max oc 5]^^^ and likewise max[9[t<;]] oc 6~^]^. Thus, as 6bl is decreased, the instability shifts to 
shorter wavelengths and becomes more rapid, while at the same time, modes having kydsL ^ 1 
are stabilized for a given value of ky. Since any real shear layer is likely to have a nonzero width, 
one may therefore remark that taking the vortex sheet limit for M ^ Merit masks the instability. 



5.2. Numerical Verification 



To independently verify our dispersion relation (C20), we ran direct hydrodynamical simula- 
tions of shear layers using the Godunov code Athena ( Stone et al.||2008 ) and compared the growth 
rate of the fastest growing mode predicted by our dispersion relation to that obtained in the sim- 



ulations. We initialized a flow along the y-direction using the setup described by Equations (75), 
and for all of the simulations we used 7 = 5/3, 5bl = 1 for the half width of the shear layer, 
and periodic boundary conditions in the y-direction; Table [T] summarizes the simulation specific 
parameters. To seed the instability, we initialized random perturbations to of magnitude 10~^ 
in the region x > —\. We found that for the M = 5 runs, we needed a high resolution in the 
x-direction to get converged estimates for the growth rates. Figure [4] shows averaged over 

box as a function of time. The solid lines show simulation results, and the dashed lines show the 
predictions from considering the fastest growing mode for the upper branch, using the dispersion 
relation (C20). There is good agreement between the two, especially for the limiting cases of e = 
and e = 1. 



Label 


M 


e 


j;-range 


y-range 


X Ny 


(BC-xl, BC-2;2) 


A 


5 


1 


(-4,4) 


(-16,16) 


8192 X 8192 


(outflow, outflow) 


B 


5 


.25 


(-4,4) 


(-16,16) 


8192 X 8192 


(outflow, outflow) 


C 


5 


.1 


(-4,4) 


(-16,16) 


8192 X 8192 


(outflow, outflow) 


D 


5 





(-1,4) 


(-16,16) 


2048 X 2048 


(reflecting, outflow) 


E 


.05 


1 


(-4,4) 


(-8,8) 


512 X 512 


(reflecting, reflecting) 



Table 1: Parameters for the simulations with Athena. M ~ Mach number, e - density ratio, x, y- 
range - box size in x, y-direction, Nx x Ny - number of cells in x, y-direction, (BC-xl, BC-x2) - 
lower and upper boundary conditions in x direction. 
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Fig. 3. — Dispersion relations for M = 5 and a finite width shear layer. The top row is for e = 1, 
the second row for e = .25, the third row for e = .1, the fourth row for e = .01, and the bottom 
row for e = 0. The first column shows the real part of ip, the second the imaginary part of ip and 
the third the imaginary part of oj. The solid line corresponds to the upper branch and the dashed 
line to the lower branch. In the top row, the upper and lower branches have the same growth rate. 
The arrows in panel (o) show the wavenumbers corresponding to trapped wavemodes (^5.3[). 
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t 



(a) 

Fig. 4. — Comparison of the analytically derived growth rates (given by the slopes of the dashed 
lines) to the ones obtained using Athena (solid lines). The curves A, B, C, D correspond to e = 1, 
.25, .1, and 0, respectively (see text for further simulation details) and have been offset vertically 
for clarity. 



5.3. Physical Intuition 

Panels, (a), (b), and (c) of Figure [s] show the spatial structure of Vx for simulations A, E &: 
D during the linear stage of the instability. The fastest growing modes in these three cases are 
very different, illustrating the different instability mechanisms which operate for the M ^> Merit 
vs. M <C Merit cases and also for e = 1 vs. e = 0. The case M <C Merit, in panel (b) has been 
extensively studied (e.g. ( Vallis|[2"006 )), so we will not consider it here.. 



For the M ^ Merit case with e = 1, the instability is caused by a radiation mechanism. This 



mechanism has already been discussed by Glatzel (1988), so we only mention it here briefly. Each 



mode of the dispersion relation can be associated with a pseudo-energy that is conserved. If a 
mode has a negative pseudo-energy, then radiation of energy away from the boundary layer region 
will cause the pseudo-energy to become even more negative, amplifying the mode and leading to 
instability. The radiation mechanism is responsible for the smooth, broad hump in panels (b) and 
(c) of Figure [3j 

For the case of M ^ Merit and e = 0, the instability mechanism is quite different. Rather than 
a broad hump, panels (n) and (o) of Figure [s] exhibit sharp localized peaks. The physical cause of 
these peaks is due to over-reflection of modes that become trapped between the lower boundary 
and the critical layer. Thus, we shall call this the over-reflection mechanism. 



The over- reflection mechanism is discussed in Narayan et al. (1987), who considered a rotating 
shear layer adjacent to a reflecting wall. They explained the instability in the e = case in terms 
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of the leaking of action density past the corotation radius, which is the analog of the critical layer 



for a rotating system. Like the pseudo-energy of Glatzel (1988), the action density is a conserved 



quantity, and instability occurs because a wavemode that is trapped between the critical layer and 
the wall has a negative action density, whereas positive action leaks out past the critical layer. As a 
result, the amplitude of the trapped wavemode becomes even more negative and grows with time. 
Put in this way, it is clear that the over-reflection and radiation mechanisms are related. However, 
they lead to a very different structure for the dispersion relation, as evidenced by panels (n) and 
(o) vs. panels (b) and (c) of Figure |3j Thus, we prefer to regard them as separate mechanisms, 
but point out that in both cases, instability is ultimately caused by radiation emitted from the 
boundary layer region. 

We now elucidate some of the properties of the modes that are trapped between the critical 
layer and the lower reflecting boundary using a simple model. Consider a sonic mode that is trapped 
between the lower boundary and the critical layer. If the mode has a well-defined phase velocity, 
then by performing a velocity boost one can transform into a frame in which the wavefronts are 
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stationary. Let us work in this frame, since it makes the explanations more clear. In order for the 
wavefronts to be stationary, the equation for a wavefront is 

dy 



±^/M{xy-l. (77) 

Here, M{x) is the Mach number as a function of x, and the positive sign is for waves propagating 
in the —x direction, while the negative sign is for waves propagating in the +x direction. Equation 



(77) is obtained from the fact that the wavefront propagates at the speed of sound, and in order for 



it to be stationary, the angle that the wavefront forms with the y-direction obeys sin(^) = 1/M(x) 



Consider now the upper branch from ^3.1, for which the critical layer is at Xc = 5bl{M — 1) /M . 
At the critical layer, we set dy/dx = 0, which is just to say that an upward traveling sound wave 
is reflected there. It immediately follows that inside the region of shear 



M{x) = -M[l-—y-5BL<x<5BL. (78) 

(79) 

Thus, the frame in which sound waves are stationary and reflect at the critical layer is boosted by 



-M relative to the frame we defined in Equations (75). 



We now consider the fate of a wavepacket that is trapped between the lower boundary and 
the critical layer. Figure [g] shows a set of wavefronts derived by integrating Equation (77) with 



red segments corresponding to propagation in the +x direction and blue segments corresponding 
to propagation in the —x direction. Consider a localized wavepacket at point A in Figure [6^. The 
wavevector of the wavepacket is oriented perpendicular to the wavefront, and as the wavepacket 
propagates towards point B, its wavevector is rotated by the shear. When the wavepacket reaches 
the critical layer at point C, it is reflected back towards point D, and the reflected wavepacket 



has a higher amplitude than the incident one (see Narayan et al. (1987) for the case of a rotating 
shear layer). After reflecting off the critical layer, the wavepacket propagates downward to point 
D, its wavevector continuing to be rotated in the same sense as before due to the shear. Finally 
the wavevector reaches the perfectly reflecting boundary at point E, whereupon the x-component 
of its wavevector is reflected and the cycle begins anew. It is thus clear that the repeating cycle 
A-E will lead to exponential amplification of the wavepacket, due to over-reflection at the critical 
layer. 

The distance in y between points A and E is the wavelength of the longest wavelength mode, 
Amax,y, that can be trapped between the lower boundary and the critical layer. A trapped mode can 
have a wavelength shorter than \me,x,y, as long as its wavelength satisfies the relation \y = Xninx,y/n 
where n is an integer. Thus, the relation for the wavenumber of the n-th trapped mode in the general 
case is 

27rn 

kn,y = T . (80) 

^max,y 
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^0 5 10 15 20 25 30 35 ^0 5 10 15 20 25 30 35 

Y Y 

(a) (b) 

Fig. 6. — The n = 1 and n = 2 trapped modes for M = 5 and 6bl = 1- The lower black line 
denotes the solid lower boundary, the upper black line denotes the top of the shear layer, and the 
dashed black line denotes the critical layer. The red segments correspond to propagation in the +x 
direction and the blue segments correspond to propagation in the —x direction. 



We indicate using arrows the values of kn,y for the first six trapped modes in panel (o) of Figure |3j 
and it is clear that they agree well with the locations of the peaks in the dispersion relation, which 
lends support for the over- reflection argument. We also plot the n = 2 case in Figure [HJd and point 
out the similarity between the analytical curves in Figure [6] and the shape of the wavefronts from 
simulation D in Figure [Sjj. 

We have discussed the radiation mechanism for e = 1 and the over-reflection mechanism for 
e = 0. We now discuss what happens in the more general case of < e < 1. As we can see from 
Figure [3j both mechanisms operate simultaneously for < e < 1. In this case, there is partial 
reflection at the lower boundary, and as e goes from one to zero, less and less of the radiation can 
escape from the boundary layer region, so the radiation mechanism becomes weaker and weaker. 
This is evidenced by the decreasing size of the bump in the second and third columns of Figure 
[3] as e goes to zero. For e = 0, the bump disappears entirely, and the radiation mechanism no 
longer operates. On the other hand, as e goes from one to zero the reflection mechanism becomes 
stronger and stronger, since more and more of the energy is reflected at the lower boundary with 
total reflection at e = 0. Interestingly, there are still some small-scale wiggles in the dispersion 
relation, even for e = 1. This is because even if the density is everywhere uniform, radiation 



can still partially reflect off the discontinuity in the velocity derivative at a; = — 1 (Glatzel 1988). 
However, reflections off a discontinuity in the velocity derivative are quite weak, so the wiggles are 
small, and the radiation mechanism is dominant. 
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6. Discussion and Conclusions 

We have studied supersonic shear instabihties that could drive the turbulence in the BLs of 
stars for which the disk is undisrupted by a magnetic field. Our study is aimed mainly at identifying 
the instabilities that lead to the formation of the BL when the disk just touches the surface of the 
star. The main result of our work is the identification of two types of instabilities that could operate 
in the BLs of such systems and had not been previously discussed in this context. 

The first is an instability of a vortex sheet at high Mach number caused by gravity. Although 
the vortex sheet is stable above a critical Mach number, the addition of a small amount of gravity 
destabilizes it. We have found that the eigenfrequencies of the dispersion relation in the limit G — )• 
acquire a purely imaginary term, which is proportional to the small parameter G. This has the 
effect of making the upper branch unstable. We now consider whether the instability of the upper 
branch in the limit G — )• is likely to be relevant during the initiation of the BL, when one may 



expect the vortex sheet approximation to be valid. Redimensionalizing Equation (74), we obtain 

uji ~ VtKe^''^i, (81) 

where Q.k is the Keplerian velocity at the surface of the star. During the initiation of the BL, 
we expect disk material to be less dense than stellar material, which means that e ^ 1. It then 



follows from Equation (81) that the characteristic growth rate of the instability is < and is 



independent of the wavenumber. 

If the BL is radially thin, as would be expected during the initiation phase, then the growth 



rate given by Equation (81 ) is small compared to the shear rate S ~ QkR*/Sbl S> ^k- Considering 
that S is the characteristic growth rate for shear instabilities, if there are other mechanisms for 
instability with growth rates proportional to S, they will quickly become dominant. 

In particular, we have demonstrated that the growth rate of the sonic instability of a finite width 
shear layer at high Mach number is proportional to S. The sonic instability is similar in nature 
to the Papaloizou-Pringle instability in that both are global instabilities and cannot be derived 
from a local analysis. There are two destabilizing mechanisms for the sonic instability. The first 
corresponds to emission of radiation and gives instability over a broad range of wavenumbers. The 
second corresponds to over-reflection of trapped modes and results in sharply peaked, disconnected 
regions of instability in fc-space. Because sonic instabilities operate on a much faster timescale than 
the gravity mechanism for a vortex sheet, this makes them an appealing candidate for the initial 
stages of boundary layer formation, when one might expect large shears to be present. 



We mention that Alexakis et al. (2004) have considered the Miles instability (Miles 1957a) in 
the context of boundary layer formation. Specifically, they invoked the Miles instability to generate 
mixing of WD and stellar material and explain the enrichment in heavy elements observed in nova 
explosions. The Miles instability was proposed as a mechanism for generating waves over water 
at low wind speeds and operates through a resonant interaction between the wind and the water 
wave. Due to this interaction, a component of the pressure perturbation is created that is in phase 



with the velocity of the air- water interface, and much hke pumping a swing, swings up the interface 
to large amplitudes. However, the Miles instability was initially introduced as a way of explaining 
the formation of waves on water at weak wind speeds, for which the air-water interface is stable 
to the KH instability. During boundary layer formation, however, large shears are generated, so 
we are not in the weak wind regime, and the Miles instability is likely to be swamped by sonic 
instabilities. 

There are two main astrophysical implications of our findings. First, we demonstrate that the 
initiation of the BL is likely to take only very short amount time after the first material from the 
disk arrives to the stellar surface. This is because the sonic instabilities that we explored in this 
paper have an extremely short growth rate, which is very weakly dependent on the density contrast 
between the disk and stellar material]^ Thus, mixing of the two fluids starts almost immediately 
after they come into contact. 

Second, given the efficiency with which the sonic instability operates, it is likely that it may 
play important role also in more developed phases of the BL evolution. As long as the BL pos- 
sesses some effective "boundaries" (e.g. sharp changes in the velocity of density behavior) and the 
gas fiow within it is supersonic the purely hydrodynamic sonic instabilities are going to operate 
in it potentially providing means for continuing mixing and angular momentum transport inside 
the layer. Future numerical calculations capable of following the nonlinear development of sonic 
instabilities should be able to address this issue. 

We are grateful to Jim Stone for useful discussions. The financial support for this work is 
provided by the Sloan Foundation and NASA grant NNX08AH87G. 
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A. Reduction of Equations 



Defining S = wdVt/dw to be the shear rate, and assuming 6bl ^ R* and fi* <^ we have 
that 2B ^ S ^ ^kR*/Sbl ^ ^K- We next note that |a;p > (Q'fa)])^ with exact equality at the 
location of a critical layer, where 5R[a)] = 0. According to Chimonas (1970), the growth rate of the 



fastest growing mode in a plane parallel stratified shear flow is of the order 



~ -s"^ 

^ 4 



(Al) 



where is the Brunt- Vaisala frequency (Equation (56)). If we assume that shear instabilities 
provide the turbulence in the BL, then we can estimate that [(Dp > 5^^. Using this estimate on 
the left hand side of Equation (15), we find that the first term, a)^, dominates the third term, 
~ 2Q.S by a factor of R^/6bl S> 1, so we can ignore the third term. Next, we compare 
the second to last term ojg/s^ to the last term 2Q.m/w on the right hand side of Equation (15). 
Defining = s^/g, dividing the second to last term by the last term, and ignoring constants 
of order unity we have CoR^/hsmQ > SR^/hgrnQ. ~ Ri^/hsSBLm- Next, we make the additional 
reasonable assumption that the fastest growing mode has mjR^ ~ ^7^) since hg is of order the 
scale height and sets the natural length scale in the problem. Continuing our line of reasoning, we 
then have R^/hs5BLfn ~ R*/Sbl S> 1. This means that the last term in equation (15) is negligible 
compared to the second to last term and can be ignored. Finally, we compare the first term, Cloj, 



and the second term, 2Bm/w, on the right hand side of Equation (14). Dividing the first term by 
the second term, using 2B ^ S, |wp > 5^ and m/R^ ~ hj^, and noting that Cl ~ hj^, we have 



This preprint was prepared with the AAS lAT^rjX macros v5.2. 
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CLUjR-t/mS > hgS/hsS ~ 1. Thus, both terms can potentially be of comparable magnitude and 
both must be retained. 



Keeping only the dominant terms for Sbl ^ R*^ using w ^ R^,, and assuming that p' j p':^ R 
and 6u'/du S> R^^, Equations (14) and (15) finally reduce to Equations (16) and (17). 



-1 



B. Generalized Rayleigh Equation and Boundary Conditions 

B.l. Derivation of the Generalized Rayleigh Equation 

We start with Equations (19) and (20). Using the definitions from ^ Equation (19) can be 
written as 



6u' + Ci{x)6u = C2{x) 

C1{X) = -{kg + W'/W) 

C2{x) = i{l-W^/s'^)ky6P/pW, 



and similarly Equation (20) can be written as 

6P' + C3{x)6P = C4(x) 



Ca{x) 



MklW^ + g{kg + ks))5u/kyW. 



(Bl) 
(B2) 
(B3) 

(B4) 
(B5) 
(B6) 



We can use the integrating factor method of solving first order differential equations to make a 
change of variables and get rid of the 6u term in Equation (Bl ) and the 5P term in Equation (B4). 
Making these changes of variable, we have 

(1 - W^/s^)ky5P 



dq 
6P' 



pW^f 



ifpjklW^ +g{kg + ks))5u 
kyW 



= {iWf)-^5u 

, 5P = f5P 



We can combine Equations (B7) and (B8) into a single equation for 

pW'^6q' 



1 - wys^ 



p{klW^ +g{kg + ks))5q. 



(B7) 
(B8) 

(B9) 



This is the same basic form as [Chimonas (1970). Since Equation (B9) is a second order differential 
equation for we can convert it to standard form (get rid of the first order term) with the change 
of variable 



W5q/ky. 



(BIO) 



Making this change of variable yields Equation (21) which is the generalized Rayleigh equation 
(Alexakis et al.||2002h. 
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B.2. Derivation of the Boundary Conditions at the Interface 

We present here a physically motivated derivation of the interfacial boundary conditions, Equa- 



tions (31) and (32). The two classical boundary conditions for the vortex sheet that should be 



satisfied at the interface are (1) that the two fluids should stay in contact at the interface, and 
(2) that the force exerted on the lower fluid by the upper fluid is equal and opposite to the force 
exerted on the upper fluid by the lower fluid. These can be expressed as 



(Bll) 
(B12) 



where 5.^ is the displacement of the interface in the x-direction, and A denotes the Lagrangian dif- 
ferential. However, we must formulate both of these boundary conditions in terms of the generalized 
streamfunction perturbation 6(j). 



We begin with the condition (5^+ = (5^- . We start with the relation 

D6^ 



Dt 



-i{uj — kyVy)5S^ = 6u, 



(B13) 



where D/Dt, denotes the Lagrangian derivative. From Equations (Bll) and (B13) we have 

UJ — kyVy bJ + kyVy 



(B14) 



Using the definitions of 5(j), W, and p given in ^ and noting that p = p at the interface (x = 0), 



we immediately get Equation (31): 



W- 



(B15) 



For the second boundary condition, we begin by writing 

AP = 6P- gp6^. 



Using Equation (B7) we can substitute for 6P in terms of 6q, which gives at the interface 



Ky 



(B16) 



(B17) 



at Equation (32): 



Using Equations ( |B13[ ), ( |B7| ), and ( |B10[ ) to substitute for 5^, 6q, and 6u in terms of 6(j), we arrive 

W-dd)'^ - Wi5(j)^ - (B18) 
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C. Finite Width Shear Layer 

We derive here the dispersion relation for a finite width shear layer in the absence of gravity 
and for a constant shear. Our analysis extends the work of Glatzel (1988) to arbitrary density 
ratios above and below the shear layer. We will use his notation in this section and consider the 
pressure perturbation 5P rather than the generalized stream function 5(f). Where appropriate, we 
describe how to transform the results back into the notation used in the body of the text. 



C.l. Setup of the Problem 



Consider the velocity profile 



and the density profile 



V{x) 



-1, X < -1 

X, —1 < X < 1 
1, X > 1, 



(CI) 



p{x) 




(C2) 



The adiabatic index and equilibrium pressure are everywhere constant so we have p_s^ = p+s^ = 
Po-Sq- The perturbations are assumed to be of the form 5P = 6 f (x) ex.p[i{kyy — ut)]. We define 
M = 1/so to be the Mach number at x = 1 inside the shear layer. We also define 



Po/p- 
Po/p+ 



(C3) 
(C4) 



and adopt the following quantities from Glatzell ( 1988 ) 



a = 






ky 


Q = 




C = 


ikyMa"^ 




1 ky 


X = 


4M 




3 


p- = 


4' 



(C5) 

(C6) 
(C7) 

(C8) 
(C9) 



Glatzel (1988) has shown that Q satisfies Whittaker's equation 



+ 



1 I X , 1/4- 

■4 + c 



p 



Q = 



(CIO) 
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and inside the shear layer has the sokition 

Q = ciM^,^(C)+C2M^,_^(C), 



(Cll) 



where ci and ci are constants and M^^^(C) is a Whittaker function. Outside the shear layer the 
velocity is constant, and the perturbation can be written as 



Q oc 



exp {±.ky(\ - eZ^M'^a'^yl'^x] , x < -I 



(C12) 



exp [±ky{l - e-^M'^a^y/'^x] , x>l 
The signs in the exponentials should be chosen based on the condition of outgoing waves (^3.1[). 



C.2. Dispersion Relation 

To derive the dispersion relation, we must apply the contact and pressure continuity boundary 



conditions at x = ±1 (Appendix B.2). Since the velocity and hence a are everywhere continuous. 



the pressure continuity boundary conditions at j; = ±1 are simply 

Qin = Qout (C13) 
where the subscripts "in" and "out" denote evaluation of a quantity inside or outside the shear 



layer, respectively. Using the expression (Glatzel 1988) 



du 



kpa 



6P' 



the contact boundary conditions at x = ±1 can be written as 



+ 



6Uin 
Pout^P'n 
1 dQin 



P'm^Pout 

_^ Pm/Pout f Pout _ 

2 V Pin C 



1/2 



(C14) 



(C15) 
(C16) 

(C17) 



4C Qin dC -z \pi 

In deriving the expression (C17), we have made use of relations ( C5||C8 ), and also the results (C12) 



and (C13). The sign in Equation (C17) should be chosen on the basis of outgoing waves, and 
the lower sign should be chosen at x = 1 and the upper sign at x = —1. We also point out that 
Pin/ Pout = e± at X = ±1. 



Next, we substitute Equation (Cll) into Equation (C17) at x = ±1, and set ci = 1, which 



sets the normalization. Using the property of the Whittaker function ( Abramowitz &: Stegun|1972 ) 
that 



c 



dC 



(C18) 
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we can derive the dispersion relation 



1 - 4x + 2C+ + 2e+C+ e 



,-1 



4JL 



1/2 



M^3(C+) + (4x + 5)M 3(C4 



1 - 4x + 2C+ + 2e+C+ e 



.-1 



4^ 



1/2 



M^^_3(C+) + (4x-l)M^^i^_3(C+) 



1 - 4x + 2C- - 2e_C- (el^ - 4^ 




M^_3(C-) + (4x + 5)M^+,^3(C-) 


l-4x + 2C--2e_C- (e='-4^) 


1/2- 


M,,_3(C-) + (4x-l)M^^,,_3(C-) 



(C19) 



Using the relation between the confluent hypergeometric function and the Whittaker function 



(Glatzel 1988), Equation (C19) can be written in terms of the confluent hypergeometric function 



as 





l-4x + 2C+ + 2e+C+(e;'-4^) 


1/2- 


i^i(i - X, i, C+) + (4x + 5)iFi(| - X, §, C+) 


l-4x + 2C+ + 2e+C+(e;'-4^)'^' 


ii"i(-i - X, -hc+) + (4x - i)ii"i(-| - X, -h c+) 





l-4x + 2C--2e_C-(el'-4 


i) 


1/2- 


iFi(|-x,i,C-) + (4x + 5)iFi(4-x,i,C-) 


l-4x + 2C--2e_C-(e:^-4^) 


1/2- 


ii^i(-i - X, -hc-) + (4x - i)ii"i(-| - X, -hc-) 



(C20) 



Equation (C20) is a generalization of the dispersion relations considered by Glatzel (1988). 
For instance, taking e-t — )• we recover his dispersion relation (5.23) for the reflecting boundary 
condition 5u{x = ±1) = 0, taking e± — t- cxd we recover his dispersion relation (5.24) for a vacuum 
boundary condition 6P{x = ±1) = 0, and taking e± = 1 we recover his dispersion relation (5.38) 
for an infinite fluid with uniform density everywhere. 



C.3. The Vortex Sheet Limit 



We can show that Equation (C20) reduces to the dispersion relation for the vortex sheet, 



(Equation (40)), by taking the limit A;^^ — t- 0. In this limit, (" ~^ and x ~^ 0, but a, the phase 
speed in the comoving frame, is finite, and 

aX 



"C 



(C21) 



is also finite. We also note the property of the hypergeometric function that in the limit C ~^ 



iFi(/,m,C) = l + -C + 0(C). 



(C22) 
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Using these relations, Equation (C20) can be reduced to the form 



(C23) 



Equation (40) can then be obtained by setting = 1, e_ = e, and using the relation M = l/s+. 



